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UV-filtered fermionic Monte Carlo 
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The short-range modes of the fermionic determinant can be absorbed in the gauge action using the loop 
expansion. The coefficients of this expansion and the zeroes of the polynomial approximating the remainder can 
be optimized by a simple, practical method. When the multiboson approach is used, this optimization results in 
a faster simulation with fewer auxiliary fields. 



Dynamical fermion simulations are orders of 
magnitude slower than quenched ones. With Hy- 
brid Monte Carlo (HMC), this remains true also 
for heavy dynamical quarks: even in the quenched 
regime, HMC is 0(100) less efficient than a lo- 
cal Monte Carlo update scheme Q. This poor 
behaviour may explain why HMC computing re- 
quirements grow rather slowly as the quark mass 
is decreased. It clearly indicates room for im- 
provement, which hopefully will persist down to 
light quarks. 

The fermion determinant can be expanded in 
loops: 

det(l-KM) = e Tr Log(l-KM) = e -E, ^Tr M' (1) 

using Wilson fermion notation. The first nonzero 
term (I — 4) gives a shift A/3 of the gauge cou- 
pling. It is remarkable that the bulk of the 
fermionic effects can be reabsorbed into this sim- 
ple shift, which can be computed perturbatively 
down to rather small quark masses [Q]. Unfortu- 
nately, neither HMC nor the alternative Multi- 
Boson (MB) method makes use of this obser- 
vation. Effective loop actions (||) 5 ff cannot 
be used to guide an exact algorithm, because 
((Scxact — S c ff) 2 ) cx Volume, which causes an 
exponentially small e - Volume Metropolis accep- 
tance. Instead, I propose to make use of the iden- 
tity e _Tr x det e A = 1 to rewrite the fermion 
determinant as 

det(l - kM) = (2) 
e l^j ' x det (1 — kM) e ^ 



The number of nonzero coefficients aj and their 
values are subject to optimization. Note that the 
first four terms cause no overhead. The new oper- 
ator whose determinant is to be taken will hope- 
fully be easier to sample than the original one. In 
particular, for very heavy quarks the coefficients 
a.,- should tend to the Taylor values n 3 /j, leaving 
only tiny higher-order terms in the determinant. 

Both HMC and MB sample the determinant 
det A by building a polynomial P n (A) which ap- 
proximates A . Gains in efficiency come from a 
better approximation (smaller degree or smaller 
error). To illustrate the advantage of using 
Eq.(2), in Fig.l I show the relative error of 
the polynomial approximation for a real variable 
x G [0.01,1]. The dotted line corresponds to 
the Chebyshev approximation used in the MB 
method; the solid line is the error achievable by 
including two terms a\ and <Z2 in the "precondi- 

tioner" e ° j x , with a polynomial of the same 
degree n — 20. Fig. 2 shows the zeroes of the two 
polynomials in the complex plane. While the ze- 
roes of the Chebyshev polynomial are distributed 
uniformly over an ellipse covering the entire ap- 
proximation interval, the other zeroes are much 
more concentrated near the origin. This is the 

essence of the benefits of the "filter" e ^j a3X ; 
large eigenvalues (corresponding to UV modes of 
the Dirac operator) are effectively removed from 
the determinant, which frees the zeroes of the 
polynomial to concentrate near small eigenval- 
ues corresponding to IR modes and genuine long- 
range fermionic interactions. 




Figure 1. Error of the Chebyshev (dotted line) 
and UV-filtered (solid line) approximations for a 
polynomial of degree 20 over the interval [0.01, 1]. 
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Figure 2. Zeroes of the Chebyshev (o) and UV- 
filtercd (+) polynomials in the complex plane. 



Choosing the coefficients {a,j} and the zeroes 
{zk} of the polynomial approximation seems like 
a stiff nonlinear optimization problem, with the 
dj appearing in the exponent. Even setting a,j = 
Vj, the choice of the zeroes {zk} has been the 
object of several studies [Hp). The problem is 



how to fix the parameters of W = Ylki^ ~ K ^ ~ 
z fe l) • (1 - kM) ■ e^=o a > AP such that detW « 
1. Rewriting the determinant as an average over 
Gaussian random vectors r\ 



det- 2 W - J dV dVe 



= l e -\Wr,\ A + \v[ 



,(3) 



we see that a sufficient condition for detVF ~ 1 
is || Wrj — 77 || 2 ~ 0. This optimization problem is 
easy to solve, using the following steps: 

1. Draw one (or more) Gaussian vectors r); 

2. Assign values to the {<Zj}; compute tp = 

(1-kM) e^7=o a i M3 m 

3. Construct the polynomial Q n (M) = YYki^ ~ 
kM — Zfel) which minimizes c =|| Q n (M)ip — r] || 2 ; 

4. Compute V a = |^f"| as a by-product; if 
|| V a || > e, return to (2). 

Step 3 is a straightforward quadratic minimiza- 
tion similar to GMRES, while the loop 2-4 uses 
Newton's method to solve the nonlinear mini- 
mization in the Oj. In principle, some subtle aver- 
aging over gauge fields should also be performed; 



in practice, different equilibrium gauge fields yield 
similar results. 

This optimization method is simple and prac- 
tical, ft requires no a priori knowledge of the 
Dirac spectrum. It is applicable to all cases where 
a determinant is involved: Wilson or staggered 
fermions, SUSY, condensed matter. When ap- 
plied to the standard MB method ({a,j} = 0, Her- 
mitian or non-Hermitian Wilson Dirac operator) , 
it reveals that the usual Chebyshev approxima- 
tion is not optimal. 

In an actual MC simulation {J7 id} — > 
{[/now}, the Metropolis acceptance probability is 

min(l, ( e -lw^'rii i +w* )v)- This can be estimated 
as erfc(c || Woid?? — rj ||) with c ~ 0(1) as a by- 
product of the above optimization. In this way, 
the benefits of adding more terms to the precon- 
ditioner (especially 6-link and larger loops) can 
be assessed without having to program the actual 
Monte Carlo update. 

Tests of the dynamical behaviour of this UV- 
filtered MB algorithm have been performed for 
moderate quark masses: (3 = 5.3, k = 0.158, 
with two flavors of Wilson fermions, on an 8 4 lat- 
tice. Table I compares the efficiency of HMC, the 
non-Hermitian MB ||, and the present method. 
All the programs used even-odd preconditioning. 
The HMC program incorporates multiple step- 
size integration [Q and low-accuracy solution dur- 
ing the trajectory j8|, and uses the BiCG75 solver 
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HMC 


MB 


this work 


A/3 








0.166 


deg. of polynom. 


26/65 


20 


7 


TintiP) in V> x v 


- 16000 


- 38000 


- 3200 



Table 1 

Comparison of three exact algorithms: Hybrid 
Monte Carlo, MultiBoson, and UV-filtered Multi- 
Boson. A/3 is the shift in the gauge coupling in- 
duced by UV-filtering. In the case of HMC, the 
degree of the polynomial is the average number of 
iterations required by the BiCG7s solver, during 
(26) and at the ends (65) of each trajectory. The 
integrated autocorrelation time for the plaquette 
is measured in units of multiplications by ip. 

|9|. The number n of fields used by the non- 
Hcrmitian MB is consistent with the number of 
iterations of the HMC solver. In contrast, the 
UV-filtered version uses ~ 3 times fewer fields. 
This solves the memory bottleneck of the MB ap- 
proach and dramatically reduces the work per in- 
dependent configuration, which grows like ~ n 2 
|9j. Fig. 3 shows the autocorrelation function of 
the plaquette for HMC and the UV-filtered MB. 
For heavier quarks the superiority of the latter 
would be even greater. For lighter quarks the ad- 
vantage is less pronounced, and details of the im- 
plementation and the tuning (eg. , over-relaxation 
of the auxiliary fields) become more relevant. 
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Figure 4. Complex spectrum of (1 — k 2 M 2 ) es- 
timated from the tridiagonal matrix generated 
by the BiCG7s solver, and zeroes of the polyno- 
mial used by the UV-filtered MB. Only one of the 
seven zeroes is devoted to controlling UV modes, 
while the other six are dedicated to the IR ones. 



